#pragma once

#include <lapacke/lapacke.h>
#include <cblas/cblas.h>
#include <vector>
#include <string>
#include <filesystem>
#include <fstream>
#include <iostream>
#include "spline.h"

namespace fs = std::filesystem;

/**
 * @brief 调用lapack求解矩阵matx=b，返回解向量x
 *
 * @param mat 系数矩阵
 * @param b 右端项
 * @return vector<double> 解向量
 */
std::vector<double> solve(std::vector<std::vector<double>> mat, std::vector<double> b)
{
    int n = mat.size();
    int m = mat[0].size();
    std::vector<double> A(n * m);
    for (int i = 0; i < n; ++i)
        for (int j = 0; j < m; ++j)
            A[j * n + i] = mat[i][j];

    int lda = n;  // Leading dimension of A (lda >= max(1, n))
    int nrhs = 1; // Number of right-hand sides
    int ldb = n;  // Leading dimension of B (ldb >= max(1, n))
    lapack_int *ipiv = new lapack_int[n];
    lapack_int n_lapack = static_cast<lapack_int>(n);
    lapack_int m_lapack = static_cast<lapack_int>(m);
    lapack_int lda_lapack = static_cast<lapack_int>(lda);
    lapack_int ldb_lapack = static_cast<lapack_int>(ldb);
    auto info = LAPACKE_dgesv(LAPACK_COL_MAJOR, n_lapack, nrhs, A.data(), lda_lapack, ipiv, b.data(), ldb_lapack);
    if (info != 0)
        throw std::runtime_error("fitCurve - DGESV");
    else
        return b;
}

// 最小二乘法,线性回归计算收敛阶
double ConvergenceOrder(const std::vector<double> &x, const std::vector<double> &y)
{
    int n = x.size();
    if (n != y.size() || n < 2)
    {
        // 处理数据点不匹配或数量不足的情况
        throw std::runtime_error("Invalid input data.");
    }

    // 构建设计矩阵X
    std::vector<double> X_data(n * 2, 1.0);
    // 构建响应矩阵Y
    std::vector<double> Y_data(n);
    for (int i = 0; i < n; ++i)
    {
        X_data[i + n] = log(x[i]);
        Y_data[i] = -log(y[i]);
    }

    // 使用 LAPACKe 和 BLAS 计算线性回归
    lapack_int info;
    info = LAPACKE_dgels(LAPACK_COL_MAJOR, 'N', n, 2, 1, X_data.data(), n, Y_data.data(), n);

    if (info != 0)
    {
        // 处理 LAPACK 调用错误
        throw std::runtime_error("ApproximateSlope - dgels");
    }

    // 返回斜率
    return Y_data[1];
}

// 计算样条dim维在x的值
template <int Dim, int Order, SplineType t>
double evaluate(const Spline<Dim, Order, t> &sp, double x, int dim = 1)
{
    int N = sp.polys[0].size() + 1;
    int flag = 0;
    if (t == B || t == CardinalB)
        flag = Order;
    if (x < sp.knots[0 + flag] || x > sp.knots[N - 1 + flag])
        throw;
    int i = 1;
    while (x > sp.knots[i + flag])
    {
        ++i;
    }
    return sp.polys[dim - 1][i - 1].evaluate(x);
}

// 创建文件夹,如果文件夹不存在
void createFilepathIfNotExist(const std::string &filepath)
{
    // 尝试创建文件夹并重新尝试打开文件
    fs::path folderPath = fs::path("txtfile");
    if (!fs::exists(folderPath))
        fs::create_directory(folderPath);
}

// 令样条输出txt文件用于绘图
template <int Dim, int Order, SplineType t>
void outputFile(const Spline<Dim, Order, t> &sp, const std::string &filename, double grid = 0.01)
{
    createFilepathIfNotExist("txtfile");
    std::ofstream output_file("txtfile/" + filename); // 打开输出文件
    if (!output_file.is_open())
    {
        throw std::runtime_error("文件夹打开失败");
    }
    std::streambuf *originalCout = std::cout.rdbuf(); // 保存原始的 std::cout 缓冲区
    std::cout.rdbuf(output_file.rdbuf());             // 重定向 std::cout 到输出文件
    int N = sp.polys[0].size() + 1;
    int flag = 0;
    if (t == B || t == CardinalB)
    {
        flag = Order;
    }
    if (Dim == 1)
    {
        for (double x = sp.knots[0 + flag]; x <= sp.knots[N - 1 + flag]; x += grid)
        {
            std::cout << x << " " << evaluate(sp, x, Dim) << std::endl;
        }
    }
    else
    {
        for (double x = sp.knots[0 + flag]; x <= sp.knots[N - 1 + flag]; x += grid)
        {
            for (int i = 1; i < Dim; ++i)
            {
                std::cout << evaluate(sp, x, i) << " ";
            }
            std::cout << evaluate(sp, x, Dim) << std::endl;
        }
    }
    std::cout.rdbuf(originalCout); // 恢复原始的 std::cout 缓冲区
    output_file.close();           // 关闭文件
}

// 输出准确函数的txt绘图文件
void outputFile(double (*f)(double), const std::string &filename, double left, double right, double grid = 0.01)
{
    createFilepathIfNotExist("txtfile");
    std::ofstream output_file("txtfile/" + filename); // 打开输出文件
    if (!output_file.is_open())
    {
        throw std::runtime_error("文件夹打开失败");
    }

    std::streambuf *originalCout = std::cout.rdbuf(); // 保存原始的 std::cout 缓冲区
    std::cout.rdbuf(output_file.rdbuf());             // 重定向 std::cout 到输出文件

    for (double x = left; x <= right; x += grid)
    {
        double y = f(x);
        std::cout << x << " " << y << std::endl;
    }

    std::cout.rdbuf(originalCout); // 恢复原始的 std::cout 缓冲区
    output_file.close();           // 关闭文件
}

// 函数重载,输出准确曲线的txt绘图文件
void outputFile(double (*fx)(double), double (*fy)(double), const std::string &filename, double left, double right, double grid = 0.01)
{
    createFilepathIfNotExist("txtfile");
    std::ofstream output_file("txtfile/" + filename); // 打开输出文件
    if (!output_file.is_open())
    {
        throw std::runtime_error("文件夹打开失败");
    }

    std::streambuf *originalCout = std::cout.rdbuf(); // 保存原始的 std::cout 缓冲区
    std::cout.rdbuf(output_file.rdbuf());             // 重定向 std::cout 到输出文件

    for (double t = left; t <= right; t += grid)
    {
        std::cout << fx(t) << " " << fy(t) << std::endl;
    }

    std::cout.rdbuf(originalCout); // 恢复原始的 std::cout 缓冲区
    output_file.close();           // 关闭文件
}

// 计算 vector<double> 的 MaxNorm
double MaxNorm(const std::vector<double> &vec)
{
    if (vec.empty())
    {
        // 如果向量为空，返回一个合适的默认值或者抛出异常
        // 这里假设返回0作为默认值
        return 0.0;
    }
    double maxNorm = std::abs(vec[0]); // 初始化为第一个元素的绝对值
    // 遍历 vector<double> 中的元素，更新最大值
    for (size_t i = 1; i < vec.size(); ++i)
    {
        double absoluteValue = std::abs(vec[i]);
        if (absoluteValue > maxNorm)
        {
            maxNorm = absoluteValue;
        }
    }
    return maxNorm;
}